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Abstract 


In this paper, a method for obtaining nonlinear stiffness coefficients in modal 
coordinates for geometrically nonlinear finite-element models is developed. The 
method requires application of a finite-element program with a geometrically non- 
linear static capability. The MSC/NASTRAN code is employed for this purpose. 

The equations of motion of a MDOF system are formulated in modal coordinates. 
A set of linear eigenvectors is used to approximate the solution of the nonlinear 
problem. The random vibration problem of the MDOF nonlinear system is then 
considered. The solutions obtained by application of two different versions of a 
stochastic linearization technique are compared with linear and exact (analytical) 
solutions in terms of root-mean-square (RMS) displacements and strains for a beam 
structure. 
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Introduction 


Development of the high speed flight vehicle technology necessitate the further 
theoretical development to understand the fatigue mechanisms and to estimate the 
service life of aerospace structures subjected to intense acoustic and thermal loads. 
Efforts to extend the performance and flight envelope of high speed aerospace ve- 
hicles have resulted in structures which may behave in a geometrically nonlinear 
fashion to the imposed loads. 

Conventional (linear) prediction techniques can lead to grossly conservative de- 
signs and provide little understanding of the nonlinear behavior. A large body 
of work exists on the prediction of geometrically nonlinear dynamic response of 
structures. All methods currently in use are typically limited by their range of 
applicability or excessive computational expense. 

Methods currently in use to predict geometrically nonlinear dynamic structural 
response include perturbation, Fokker- Plank- Kolmogorov (F-P-K), Monte Carlo 
simulation and stochastic linearization techniques. Perturbation techniques are lim- 
ited to weak geometric nonlinearities. The F-P-K approach 1,2 yields exact solutions, 
but can only be applied to simple mechanical systems. Monte Carlo simulation 
is the most general method, but computational expense limits its applicability to 
rather simple structures. Finally, stochastic linearization methods (e.g. equivalent 
linearization, see 2-6 ) have seen the most broad application for prediction of geomet- 
rically nonlinear dynamic response because of their ability to accurately capture the 
response statistics over a wide range of response levels while maintaining a relatively 
light computational burden. 

The equations of motion of a MDOF, viscously damped geometrically nonlinear 
system can be written in the form 

MX + CX + KX + , (X) = F (1) 

where M, (7, K , X and F are the mass, damping, stiffness matrices, displacement 
response vector and the force excitation vector respectively. The vector function 
, ( X ) generally includes 2nd and 3rd order terms in X. There exist mathematical 
difficulties in the derivation of a general solution to equation (1) for the case of 
random excitation. An approximate solution can be achieved by formation of an 
equivalent linear system: 


MX + CX + (K + K e )X = F (2) 

where K e is the equivalent linear stiffness matrix. Assuming the Gaussian zero-mean 
loading and response, the convential stochastic linearization technique based on the 
force-error minimization (see for example Roberts et al. 3 , Atalik et al. 4 ) yields the 
following expression for the equivalent stiffness term (which replace the nonlinear 
one): 


K e = E[-M-\ 
l <9X j 


where E[...] is the operator of expectation of a random quanti 


( 3 ) 


2 



Elishakoff et al. 5,6 proposed alternative stochastic linearization approach, based 
on the potential energy error minimization and numerical results were demonstrated 
for the case of SDOF systems. The development of this approach for the case of 
MDOF systems was shown in Muravyov et al 7 . Namely, the following system of N 2 
linear equations for the equivalent stiffness matrix K e was obtained: 

N 

Y K ei jE[xiXjX k xi\ = 2 E[x k xilJ (X)] (4) 

i,j = 1 

where k, l = 1, N and the symmetry of the matrix K e is assumed. 

The matrix of the system in equation (4) involves 4th order moments of displace- 
ments and the right-hand side (assuming that the potential energy is a function of 
the 2nd, 3rd and 4th order terms) involves moments of 4th, 5th and 6th order. Using 
the Gaussian distributed, zero-mean response assumption means that the odd order 
moments are zero and the higher even order statistical moments can be expressed 
in terms of the 2nd order moments, e.g., 

E[xiXjX k xi] = E[x t x 3 ]E[x k xi] + E[x t x k \E[x 3 X[] + E[x t xi]E[x 3 x k \ (5) 

and for 6th order moments, it is 

E[x t x 3 x k xix m x n \ = E[x l x J \E[x k xi]E[x m x n ] + E[x l x k ]E[x J xi]E[x m x n ] + ...+ 

E[xiX n \E[xjX m ]E[x k xi\ 

where only three terms out of the 15 ones obtained by index permutation are shown. 
Therefore the matrix and right-hand side of (4) can be determined solely by the 
covariance matrix of displacements. 

Having the equivalent linear stiffness matrix defined through either the force 
error or potential energy error minimization techniques, one can proceed with the 
solution of the equivalent linear system. Assuming stationary excitation, a station- 
ary response is sought precluding the need for initial conditions. As the equivalent 
stiffness matrix K e is a function of the unknown displacement response vector, the 
solution to the system of equations of motion takes an iterative form. The descrip- 
tion of such a procedure can be found, e.g. in Roberts et al. 3 , Muravyov et al. 7 . 

The analytical form of the nonlinear terms facilitates the solution of equations 
(1) when the forces and displacements are random functions of time. In the section 
below a method for derivation of the analytical expression (quadratic and cubic 
polynomials) for the nonlinear term , (X) will be shown. 

Determination of Nonlinear Stiffness Coefficients 


One method of determining the nonlinear stiffness coefficients is through the 
use of a finite element approach. Existing finite element commercial programs are 
unable to provide these nonlinear stiffness coefficients directly. It is desirable to 
achieve a solution within a commercial finite element code to take advantage of the 
comprehensive element library, etc. necessary for modeling complex structures. This 
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section describes a method of determining the nonlinear stiffness coefficients through 
the use of the nonlinear static solution capability that exists in many commercial 
finite element codes. 

For MDOF structures, it is expedient to seek a solution in modal coordinate 
space 

X = $q (6) 

where $ is generally a subset (L < N) of the linear eigenvectors (normal modes). 
Such a representation allows the size of the problem to be significantly reduced 
without a noticeable loss of accuracy in many cases. 

From (1), one can obtain the following set of differential equations in terms of 
modal coordinates (i = 1,T): 


L 

W + kiQi(t) + ?2, ?l) = fi(t) (7) 

j=i 

where the nonlinear terms will be represented in the following form 

L L 

7 .'(? i >? 2 , -,? l ) = J2 a )kWk + J2 VjkiWMi ( 8 ) 

j,k=j j,k=j,i=k 

where the first index j takes values 1,2,...,L, the index k takes values from j (the 
current first index value) and up (j + 1 j + 2 ... to L ) and the third index / takes 
values from k (the current second index value) and up (k + 1, k + 2 to L). 

A procedure for determination of the coefficients a l jk and V jkl requires the appli- 
cation of a finite element program with a nonlinear static solution capability. In this 
study, the MSC/PATRAN and MSC/NASTRAN programs 8,9 are utilized. 

The suggested technique is based on the restoration of nodal applied forces from 
enforced nodal displacements prescribed to the whole structure in a static solution 
(linear and nonlinear). Namely, by prescribing the physical nodal displacements 
(vector X c ) to the structure, one can restore the nodal forces Ft and the corre- 
sponding nonlinear contribution F c : 

F c = , (X c ) = F t - I<X c (9) 

The displacements X c can be prescribed by creating a displacement constraint set for 
the model in MSC/PATRAN, then the nodal applied forces Ft will arise as single- 
point-constraint forces in a MSC/NASTRAN nonlinear static solution. In general, 
any set of displacements can be prescribed, e.g. for the bending problem of thin 
plate/beam-like structures, one can prescribe both out-of-plane and in-plane (mem- 
brane) components of displacement. The present analysis, however, is restricted to 
problems where in-plane displacements will be assumed small (thus neglected) and 
only the out-of-plane nonlinear stiffness coefficients are determined. 

To illustrate the technique, one can begin with the prescription of displacements 
for the whole structure in the following form 

x c = (Mi (10) 
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The nodal force vectors Ft (nonlinear static solution) and KX C (linear static solu- 
tion) are provided by MSC/NASTRAN. The nonlinear term F c can then be evalu- 
ated by equation (9). The vector of modal forces F c = F c is calculated and it is 
represented as 

F c = $ t F c = $ T , (X c ) = $ T , (<Mi) = [a* n ]9i9i + Wmlqiqiqi (H) 

where the sought stiffness coefficients [an], [^ 111 ] are column- vectors L X 1 (i = 1, L). 
Note that all other nonlinear terms in (11) do not appear since qj = 0 for j ^ 1. 

Prescribing a displacement held with opposite sign X c = —<f>\q\ results in a 
modal force vector (denoted by F- c ): 

F_ c = $ t F_ c = $ T , (X c ) = $ T , (-<Mi) = [a^qiqi - [b^qiqiqi (12) 

where the quadratic (even) term will be the same as in (11) and the cubic (odd) 
term takes on a sign change. 

Note that in the system of equations (11) and (12), the value of qi is given. The 
coefficients [a* n ] , [&* m ] (i = 1, L) can be determined from this system of 2 X L linear 
equations. In an analogous manner, i.e. prescribing X c = 4>jqj, all other coefficients 
\ a )j\ i [frjjj] can be determined. 

A similar technique can be employed to determine coefficients with two or three 
inequal lower indices, e.g., [«i 2 ] , [^ 112] 5 [^ 122 ] or [^ 123 ] • Note that coefficients of the 
latter type appear only if the number of retained eigenvectors L in (6) is greater than 
or equal to 3. Determination of coefficients [«i 2 ] , [^ 112 ] an( l [^ 122 ] w bl be considered 
as an example. Prescribe the displacement held to the model in the following form 

X c = <j>i qi + <(>292 

then the calculated (using MSC/NASTRAN) corresponding modal force vector F c 
is represented as follows 

F c = <$> T F c = $ T , (</i qi + ^ 292 ) = [ a n ]qiqi + Wmlqiqiqi + [F 22 \q2q2 + [^222]?292?2+ 

[a\ 2 \qiq 2 + [&i 12 ]?i?i ?2 + [ 6 i 22 ]? 2 ? 2 ?i (13) 

Prescribing the opposite sign displacement held 

x c = -<f>iqi - hq 2 

one obtains a second set of equations 

F-c = $ T F_ C = $ T , (-(Mi - <t> 2 q2 ) = [«' n ]9i9i - \b\n\qiqiqi + 

[b l 222 1 \ ( M 2 q 2 + [ai 2 ]? 1 ?2 - [ 6 i 12 ]?i?i ?2 - [&i 22 ]? 2 ? 2 ?i (14) 

Summing (13) and (14), one obtains 

F c + F-c = 2[a\ 1 ]qiqi + 2[a 22 ]q 2 q2 + 2[a\ 2 \qiq 2 
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From this equation, the coefficients [a* 12 ] are determined (note that the coefficients 
[a*n], [a 22 ] were already determined above). 

Now we have two sets of L equations (13) and (14), but to determine cubic coeffi- 
cients [b\ 12 ] and [b\ 22 ] from them is not possible since the system matrix has linearly 
dependent rows. Therefore, an additional type of displacement held is required. One 
can prescribe the following type 


X a — <f>iqi — 

Then the modal force vector is equal to 

F a = $ T F a = $ T , (faq i - <f> 2 q 2 ) = [aii]?i?i + [&m]9i9i9i + [a’ 22 ]9292 “ [& 222 ]? 2 ? 2 ? 2 - 

[a\ 2 \q 1 q 2 - [&u 2 ] 9 i 9 i 92 + [ 6 i 22 ]? 2 ? 2 ?i (15) 

From the system of 2 * L linear equations (13) and (15), the coefficients [b\ 12 ] and 
[^ 122 ] can be determined. In a similar manner, all coefficients of the type [&*• - fc ] and 
[b\ k j] can be determined. 

Now one can proceed with determination of coefficients with three different lower 
indices, like [u>) 23 ]. Prescribe a displacement held to the model in the following form 

x b = <Mi + hq2 + 

Calculated (by MSC/NASTRAN) the nonlinear nodal force vector F b is 

Fb = , (Xb) = , (<^i<7i + <f>2q2 + ^ 393 ) 
and the corresponding modal force vector will be 

F b = Fb = $ T , (<j)i qi + <fi 2 q 2 + ^ 393 ) = [ a n ]qiqi + [^m]9i9i9i + [F 22 \q 2 q 2 

+ [^2221^2^2^2 + [F 33 ]q3q3 + [^333]939393 + [«l 2 ]9l92 + [&112]9l9l92 + [&22l] 929291 + 
[«1 3 ]9 i 93 + [^113l9l9l93 + [^33ll93939l + [^ 23 ! 9293 + [^ 223 ] 929293 + [^ 332 ! 939392 + [&123]9l9293 

Note that all coefficients in this equation have been already determined, except 
[^ 123 ] • In analogous manner all coefficients of type (j ^ k,j ^ l,k ^ l) can be 
found. 

Having the modal equations of motion (7) formulated, solution to these equations 
can now be undertaken through a variety of techniques. For the case of random 
loading, the application of the equivalent stochastic linearization was implemented 
in this study. Within the framework of the force-based technique, the equivalent 
stiffness matrix (according to the formula (3)) will have the following form 

K = pr d(7l>72,-,7 l) h 
e [ d{ qi ,q 2 ,...,q L ) 1 

Note that the derivatives and expectations in (16) can be easily evaluated due to 
the analytical representation of the nonlinear terms in (8). 
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Derivation of the potential energy function coefficients 


In order to apply the energy-based version of stochastic linearization, it is neces- 
sary to have the expression of potential energy U of the system. Using expressions 
(8), one can proceed with the determination of the potential energy in terms of modal 
coordinates. It is known that nonlinear elastic force terms satisfy the following 

dll 

7.'(?i>?2, -,?l) = w— i = l,L (17) 

Oq t 

where U is the potential energy generated by nonlinear terms only. Since all nonlin- 
ear coefficients in 7;(gi, q 2} ...,7 l) have been determined, the potential energy func- 
tion U(qi } q 2} ..., qh) can be derived and it can be used in the energy-based stochastic 
linearization technique. 

Represent the potential energy (nonlinear term contribution) of the system in 
the following form 

L 

U = Y d s3k iq s q 3 q k qi (18) 

s,j,k,l 

where the first index s takes values 1,2,...,L, the index j takes values from s (the 
current first index value) and up (5 + 1,5 + 2 ... to L), the third index k takes values 
from j (the current second index value) and up (j + 1, j + 2 to L), and the forth 
index / takes values from k (the current third index value) and up (k + 1, k + 2 to 
L). The task now is to find coefficients d sjk i- Substituting (18) in (17) one obtains 

L . d L 

7 t(?i, ?2, •••, ?l) = Y fykiWkqi = ^-( Y d sjkiq s q 3 qkqi ) (19) 

j,k,l 9q i s,j,k,l 

Note that the quadratic terms in functions 7;(gi, q 2} were omitted, otherwise 

it would be necessary to introduce cubic terms in the expression (18). The quadratic 
terms in 7;(gi, q 2} ..., qi) arise in problems where membrane (in-plane) displacements 
of bending plate, beam-like structures occur. In this study the analysis is restricted 
to the flexural vibration problem with no in-plane components of the displacement. 
From (19) the sought coefficients d S]k i can be found as follows 


d 


sjkl 


' b %i/ 2 if S =J 
< b)ki / 3 iis = j = k 

Hs=j=k=l 
b s ]kl otherwise 


which can be easily verified by substitution of U(q\ } ...,5l) in the expressions (17). 

Now when the expression for U in terms of modal coordinates has been de- 
rived one can proceed with application of energy-based version of linearization (see 
definition of the equivalent stiffness matrix in (4) for this case). 
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RMS Strain Calculation 


When the solution of nonlinear equations of motion (7) is accomplished and the 
covariance matrix of displacements in physical coordinates is recovered, one can 
proceed with the determination of the variance/rms of strain components. Consider 
such a procedure of the strain determination for a flexural vibration problem of a 
beam in plane X — Z (Fig. 1). The Green strain component e x at the point with 
the coordinates (x,z) is 


e x (x,z) 


d 2 w(x) 1 
Z dx 2 + 2 


dw(x) 2 
dx 


( 20 ) 


where w is the flexural displacement and the in-plane (membrane) component of 
displacement is not taken into account, since it was assumed that (for certain prob- 
lems) it is small and can be neglected. Note that only the out-of-plane nonlinear 
stiffness coefficients were introduced in this study. 

It is necessary to use the shape functions to represent the displacement and 
strain at any location of the beam. Considering two coordinates for each node 
w and 6 (rotation about ffi-axis), one can introduce four shape functions for each 
beam element, namely, fi(x) } f 2 (x) associated with the 1st element node coordinates 
uq, 6 1 and f 3 (x) } f 4 (x) associated with the 2nd element node coordinates w 2} 6 2 . 
These shape functions fi(x) } f 2 (x),f 3 (x), f 4 (x) can be found, for example, in Weaver 
et al. 10 . Note that they are not necessary coincide with the NASTRAN’s shape 
functions of the beam type element, but since the NASTRAN’s shape functions are 
not available, these functions were employed. 

Therefore assuming that for points within each element the flexural displacement 
is 

w = wifi(x) + 0if 2 (x) + w 2 f 3 (x) + e 2 f 4 (x) 
from (20) one obtains 

e x (x,z) = -z(fl(x)w 1 + f 2 (x)0 1 + f 3 (x)w 2 + fl{x)0 2 ) + ^(f' 1 (x)w 1 + /'(x)ffi + 


f' 3 (x)w 2 + f' 4 (x)9 2 ) 2 

and the expectation of e 2 will be 

4 ] = ^(/;/;£[^+2/;/;a[^ 

2/^/sE^w,] + 2fZf';E[e 1 e 2 ] + f’’f’ 3 E[w 2 2 ] + 2f’’f:E[w 2 9 2 ] + /T/MH 
\ E [{f'l{ X ) W l + /2(®)01 + fii X ) W 2 + f4i x )0 2 ) 4 ] 

where the 3rd moments were omitted, since the Gaussian, zero-mean response is 
assumed. The terms with the co-factor | produce moments of the 4th order which 
in turn, by (5), are expressed in terms of the 2nd order moments. Thus using 
the covariance matrix of displacements and shape functions of elements, one can 
calculate the variance of strain components (in this example E[e 2 x ]) at any location 
of the structure and the root-mean-square (RMS) strain \J E[t 2 x ] can be calculated 
as well. 
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Young’s 

modulus 

Poisson’s 

ratio 

density 

0.73e+ll 

0.325 

0.2763e+04 

length 

width 

thickness 

0.4572 

0.0254 

0.002261 


Table 1: Parameters of the beam model 

Numerical Results for a MDOF beam structure 

The method which employs the technique described above, i.e. the determination 
of the nonlinear stiffness coefficients plus the two versions of stochastic linearization 
was implemented in a new in-house code. 

The numerical results presented in this section correspond to a clamped-clamped 
beam model in Figure 1. The finite-element model had 19 nodes and 18 beam 
elements. The parameters of the model in Fig. 1 are shown in Table 1 (system of 
units is SI, [to], [N/m 2 ] } [kg] etc.), where the width and thickness are dimensions 
of the cross-section of the beam. The first two natural frequencies (associated with 
flexural symmetric modes in the excitation plane) are 57.4 Hz and 310.1 Hz. These 
two mass-normalized modes were chosen to approximate the motion according to 
formula (6). 

The nonlinear stiffness coefficients determined with application of the procedure 
described above are summarized in Table 2. The quadratic terms were negligible, 
so only the 3rd order terms are shown. Since the modal coordinates q\,q 2 are 
nondimensional, the units of these nonlinear coefficients are in [N * to]. 

Note that from (17), it would follow that 

dj l _ dj k_ _ d 2 U 
dq k dqj dq k q 3 

Comparing the terms with like powers in qj and q k leads to a certain relation between 
the nonlinear coefficients, for example, for the cubic coefficients b\ 22 and b\ 12 it is 


and for other types, it is 


orl _ i2 ol2 _ 7l 

°"222 ~ "122 °"lll — "112 

It turned out that the computed nonlinear stiffness coefficients (see Table 2) are in 
an excellent agreement with these relations. 

The RMS flexural displacements and RMS strain component e x of the points 
located on the surface layer of the beam and along the coordinate x were obtained 
by using the force-based and energy-based versions of stochastic linearization and 
results are compared with the linear solution and with exact (FPK) solution for 
some cases below. 
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!>ln 

^222 


^122 

0.899e+12 

0.977e+13 

0.191e+13 

0.139e+14 

f>ni 

b 2 

"222 

blu 

b 2 

"122 

0.638e+12 

0.608e+14 

0.139e+14 

0.293e+14 


Table 2: Nonlinear stiffness coefficients for the beam model 


The system stiffness and damping matrices of the beam model in modal mass- 
normalized coordinates were as follows 


' 1. 3009877 + 05 

0 

O = 

' 4.039 

0 

0 

3.79653 E + 06 

0 

4.039 


At first, the following type of the spectral density matrix of the external forces (in 
modal coordinates) was assumed: 


Sf(co) 


a 0 

0 a 


which constituted a white noise excitation with the frequency band taken from 0 to 
550 Hz. For this type of excitation the FPK solution was possible to obtain, see, 
for example Bolotin 1 . The RMS flexural displacements as function of coordinate x 
are shown in Fig. 2, where due to the symmetry only a half of the beam is shown. 
One can see that the force-based version slightly underpredicts the values of RMS 
displacements and that the linear analysis overpredicts it by about 50 % (at the 
middle node 10). RMS displacements at the middle node for the nonlinear analyses 
achieve about 60 % of the beam’s thickness. 

The RMS strains shown in Figures 3 and 4 indicate that FPK solution is between 
the force-based and energy-based solutions for all points along the beam and for both 
cases of loading. For the linear analysis (with the term , (X)=0 in (1)) the linear 
strain recovery was applied, i.e. the 2nd term in (20) was omitted. 

As a 2nd type of excitation, a uniform (in space) white noise type pressure was 
applied to the beam. The frequency band was the same 0-550 Hz. This type of 
excitation corresponded to the following spectral density matrix of forces in modal 
coordinates: 

0.0536 -0.0236 

-0.0236 0.01052 

where the units are in [ ( JV *™Pffi* m ) ] - ppg FPK solution was not available for this type 
of loading. The RMS strains for this case are shown in Fig. 5. One can see a larger 
difference between the two versions of linearization. The linear analysis predicts 
larger values of maximum strains than the nonlinear analyses and underpredicts 
strains for some portion of the beam. The difference in maximum RMS strain (at 
the clamping) between the linear and nonlinear analyses achieves about 50 %. 

One can see that in all cases the energy-based version of linearization yields 
greater values of strain (achieving 15 % difference in some cases) than the values 


Sf(co) = 
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from the force-based version, i.e. it produces more conservative estimate for RMS 
strains. 

Conclusions 


A new method for determination of nonlinear stiffness coefficients has been sug- 
gested which utilizes a finite-element commercial software with geometrically non- 
linear static capability. It has been shown that application of MSC/NASTRAN and 
MSC/PATRAN for this purpose is sufficient. 

This method has been incorporated into a program which calculates a steady- 
state response of a MDOF structure to a Gaussian zero-mean random excitation. 
Efforts are presently underway to implement this method into MSC/NASTRAN 
through a DMAP Alter. 

Two versions of the stochastic linearization have been compared on example 
of beam structure and the results have shown that the energy-based version (in all 
cases) provided a more conservative estimate of strains than the conventional (force- 
based) version. This is a useful fact, since application of the energy-based version 
would lead to more safe estimate of the fatigue life of the structure. 
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Figure 2: RMS displacement w vs. longitudinal coordinate of the beam 
for a=0.02 ' 
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Figure 3: RMS strain t x vs. longitudinal coordinate of the beam 
for a=0.02 

raa/ s 


13 





1.0E-3 


9.0E-4 - 

8.0E-4 

7.0E-4 

6.0E-4 

5.0E-4 

4.0E-4 

3.0E-4 

2.0E-4 

1.0E-4 t 

O.OEO 



force-based 

energy-based 

FPK 

Linear 



5 6 7 

Node 


Figure 4: RMS strain t x vs. longitudinal coordinate of the beam 
for a=0.08 (JV * m) * ( / JV * m) 
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Figure 5: RMS strain e x vs. longitudinal coordinate of the beam. 
Uniform (in space) white-noise pressure 
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